WASP.version<-8

Introduction

Harpeth TMDL – Water Quality Analysis Simulation Program (WASP version 8)

Exploration of water quality data using the calibrated WASP models provided to TDEC and others by the EPA Region 4

This file was last updated on 2022-08-09 11:29:48


#Import Data


Grab <- read_excel("GrabSamples.xlsx", 
    col_types = c("skip", "text", "text", 
        "skip", "skip", "date", "text", "numeric", 
        "text"))

head(Grab,15)
## # A tibble: 15 × 6
##    DataSource Station_ID             Date_Time           Pcode      Result Units
##    <chr>      <chr>                  <dttm>              <chr>       <dbl> <chr>
##  1 STORET     11NPSWRD_WQX_NATR_BUBR 2012-01-03 08:33:00 COND     306      uS/c…
##  2 STORET     11NPSWRD_WQX_NATR_BUBR 2012-01-03 08:33:00 DOSAT     90.3    %    
##  3 STORET     11NPSWRD_WQX_NATR_BUBR 2012-01-03 08:33:00 DO        12.6    mg/l 
##  4 STORET     11NPSWRD_WQX_NATR_BUBR 2012-01-03 08:33:00 PH         7.76   std …
##  5 STORET     11NPSWRD_WQX_NATR_BUBR 2012-01-03 08:33:00 NO3_N      3.03   mg/l 
##  6 STORET     11NPSWRD_WQX_NATR_BUBR 2012-01-03 08:33:00 TEMPWC     0.4    degC 
##  7 STORET     11NPSWRD_WQX_NATR_BUBR 2012-01-03 08:33:00 FLOW_CMS   0.0212 cms  
##  8 TDEC       TNW000002791           2012-01-03 11:55:00 COND     376      uS/c…
##  9 TDEC       TNW000002790           2012-01-03 10:45:00 COND     373      uS/c…
## 10 TDEC       TNW000003300           2012-01-03 09:35:00 COND     295      uS/c…
## 11 TDEC       TNW000005654           2012-01-03 12:20:00 COND     246      uS/c…
## 12 TDEC       TNW000003508           2012-01-03 11:15:00 COND     242      uS/c…
## 13 TDEC       TNW000006704           2012-01-03 13:15:00 COND     234      uS/c…
## 14 TDEC       TNW000006704           2012-01-03 13:15:00 DO        16.5    mg/l 
## 15 TDEC       TNW000003508           2012-01-03 11:15:00 DO        15.6    mg/l

BaseModel <- read_delim("Epa_Base_Model_2011_To_2020.txt", 
     delim = "\t", escape_double = FALSE, 
     col_types = cols(`Date-Time` = col_datetime(format = "%m/%d/%Y %H:%M"), 
         CBOD = col_number(), DO = col_number(), 
         FLOW_CMS = col_number(), NH3_N = col_number(), 
         NO3O2 = col_number(), PHYTO = col_number(), 
         SOD_T = col_number(), `SOLID-WS` = col_number(), 
         TN = col_number(), TP = col_number(), 
         WTEMP = col_number()), trim_ws = TRUE)

head(BaseModel,15)
## # A tibble: 15 × 13
##    `Date-Time`         Station_ID   CBOD    DO FLOW_CMS  NH3_N NO3O2 PHYTO SOD_T
##    <dttm>              <chr>       <dbl> <dbl>    <dbl>  <dbl> <dbl> <dbl> <dbl>
##  1 2011-01-01 00:00:00 ARKANSAS_1… 0      0       0.202 0      0         0  0   
##  2 2011-01-02 00:04:00 ARKANSAS_1… 0.742  7.22    0.200 0.0116 0.153     0  1.22
##  3 2011-01-03 00:01:00 ARKANSAS_1… 0.742 10.7     0.200 0.0223 0.158     0  1.26
##  4 2011-01-04 00:04:00 ARKANSAS_1… 0.742 10.7     0.200 0.0354 0.160     0  1.25
##  5 2011-01-05 00:08:00 ARKANSAS_1… 0.742 10.7     0.200 0.0465 0.161     0  1.24
##  6 2011-01-06 00:01:00 ARKANSAS_1… 0.742 10.7     0.200 0.0550 0.161     0  1.22
##  7 2011-01-07 00:04:00 ARKANSAS_1… 0.742 10.7     0.200 0.0612 0.161     0  1.21
##  8 2011-01-08 00:08:00 ARKANSAS_1… 0.742 10.7     0.200 0.0658 0.161     0  1.20
##  9 2011-01-09 00:01:00 ARKANSAS_1… 0.742 10.7     0.200 0.0690 0.161     0  1.19
## 10 2011-01-10 00:04:00 ARKANSAS_1… 0.742 10.7     0.200 0.0711 0.161     0  1.18
## 11 2011-01-11 00:07:00 ARKANSAS_1… 0.742 10.7     0.200 0.0724 0.161     0  1.18
## 12 2011-01-12 00:01:00 ARKANSAS_1… 0.742 10.7     0.200 0.0731 0.161     0  1.17
## 13 2011-01-13 00:04:00 ARKANSAS_1… 0.742 10.7     0.200 0.0733 0.161     0  1.16
## 14 2011-01-14 00:07:00 ARKANSAS_1… 0.742 10.7     0.200 0.0731 0.161     0  1.15
## 15 2011-01-15 00:00:00 ARKANSAS_1… 0.742 10.7     0.200 0.0727 0.161     0  1.15
## # … with 4 more variables: `SOLID-WS` <dbl>, TN <dbl>, TP <dbl>, WTEMP <dbl>



Cont <- read_csv("Cont.csv", col_types = cols(...1 = col_skip(), 
    Date_Time = col_datetime(format = "%Y-%m-%d %H:%M:%S"), 
    Result = col_number()))

head(Cont,15)
## # A tibble: 15 × 5
##    Station_ID    Date_Time           Pcode  Result Units
##    <chr>         <dttm>              <chr>   <dbl> <chr>
##  1 USGS_03432390 2012-01-26 00:00:00 TEMPWC   11.8 degC 
##  2 USGS_03432390 2012-01-26 00:00:00 COND    550   <NA> 
##  3 USGS_03432390 2012-01-26 06:00:00 TEMPWC   12.6 degC 
##  4 USGS_03432390 2012-01-26 06:00:00 COND    552   <NA> 
##  5 USGS_03432390 2012-01-26 12:00:00 TEMPWC   13.6 degC 
##  6 USGS_03432390 2012-01-26 12:00:00 COND    499   <NA> 
##  7 USGS_03432390 2012-01-26 18:00:00 TEMPWC   14   degC 
##  8 USGS_03432390 2012-01-26 18:00:00 COND    451   <NA> 
##  9 USGS_03432390 2012-01-27 00:00:00 TEMPWC   12.5 degC 
## 10 USGS_03432390 2012-01-27 00:00:00 COND    306   <NA> 
## 11 USGS_03432390 2012-01-27 06:00:00 TEMPWC   11.9 degC 
## 12 USGS_03432390 2012-01-27 06:00:00 COND    455   <NA> 
## 13 USGS_03432390 2012-01-27 12:00:00 TEMPWC   12.2 degC 
## 14 USGS_03432390 2012-01-27 12:00:00 COND    512   <NA> 
## 15 USGS_03432390 2012-01-27 18:00:00 TEMPWC   12.6 degC



Natural <- read_delim("NATURAL_MODEL_2011_To_2020.txt", 
    delim = "\t", escape_double = FALSE, 
    col_types = cols(`Date-Time` = col_datetime(format = "%m/%d/%Y %H:%M"), 
        CBOD90 = col_number(), DIP = col_number(), 
        DO = col_number(), DON = col_number(), 
        DOP = col_number(), FLOW_CMS = col_number(), 
        NH3_N = col_number(), NO3O2 = col_number(), 
        PHYTO = col_number(), PON = col_number(), 
        SOD_T = col_number(), TKN = col_number(), 
        TN = col_number(), TP = col_number(), 
        WTEMP = col_number()), trim_ws = TRUE)

head(Natural,15)
## # A tibble: 15 × 17
##    `Date-Time`         Station_ID   CBOD90    DIP    DO    DON    DOP FLOW_CMS
##    <dttm>              <chr>         <dbl>  <dbl> <dbl>  <dbl>  <dbl>    <dbl>
##  1 2011-01-01 00:00:00 ARKANSAS_131   0    0       0    0      0         0.210
##  2 2011-01-02 00:06:00 ARKANSAS_131   2.38 0.0968  7.18 0.0193 0.0107    0.208
##  3 2011-01-03 00:04:00 ARKANSAS_131   2.38 0.0948 10.7  0.0193 0.0107    0.208
##  4 2011-01-04 00:02:00 ARKANSAS_131   2.38 0.0949 10.6  0.0193 0.0107    0.208
##  5 2011-01-04 23:59:00 ARKANSAS_131   2.38 0.0951 10.6  0.0193 0.0107    0.208
##  6 2011-01-06 00:05:00 ARKANSAS_131   2.38 0.0953 10.6  0.0193 0.0107    0.208
##  7 2011-01-07 00:02:00 ARKANSAS_131   2.38 0.0954 10.6  0.0193 0.0107    0.208
##  8 2011-01-08 00:00:00 ARKANSAS_131   2.38 0.0956 10.6  0.0193 0.0107    0.208
##  9 2011-01-09 00:06:00 ARKANSAS_131   2.38 0.0957 10.6  0.0193 0.0107    0.208
## 10 2011-01-10 00:03:00 ARKANSAS_131   2.38 0.0958 10.6  0.0193 0.0107    0.208
## 11 2011-01-11 00:00:00 ARKANSAS_131   2.38 0.0960 10.6  0.0193 0.0107    0.208
## 12 2011-01-12 00:06:00 ARKANSAS_131   2.38 0.0961 10.6  0.0193 0.0107    0.208
## 13 2011-01-13 00:04:00 ARKANSAS_131   2.38 0.0962 10.6  0.0193 0.0107    0.208
## 14 2011-01-14 00:01:00 ARKANSAS_131   2.38 0.0963 10.6  0.0193 0.0107    0.208
## 15 2011-01-14 23:58:00 ARKANSAS_131   2.38 0.0964 10.6  0.0193 0.0107    0.208
## # … with 9 more variables: NH3_N <dbl>, NO3O2 <dbl>, PHYTO <dbl>, PON <dbl>,
## #   SOD_T <dbl>, TKN <dbl>, TN <dbl>, TP <dbl>, WTEMP <dbl>


Base_PS <- read_delim("Harpeth_Epa_Base_Model_2_25_X_Ps.txt", 
    delim = "\t", escape_double = FALSE, 
    col_types = cols(`Date-Time` = col_datetime(format = "%m/%d/%Y %H:%M"), 
        Station_ID = col_character(), CBOD = col_number(), 
        DIP = col_number(), DO = col_number(), 
        DON = col_number(), DOP = col_number(), 
        FLOW_CMS = col_number(), NH3_N = col_number(), 
        NO3O2 = col_number(), PHYTO = col_number(), 
        PON = col_number(), SOD_T = col_number(), 
        TN = col_number(), TP = col_number(), 
        WTEMP = col_number()), trim_ws = TRUE)
head(Base_PS)
## # A tibble: 6 × 16
##   `Date-Time`         Station_ID  CBOD   DIP    DO    DON    DOP FLOW_CMS  NH3_N
##   <dttm>              <chr>      <dbl> <dbl> <dbl>  <dbl>  <dbl>    <dbl>  <dbl>
## 1 2013-01-01 00:00:00 ARKANSAS_…  0    0      0    0      0        0.0421 0     
## 2 2013-01-01 06:03:00 ARKANSAS_…  1.87 0.431  5.07 0.114  0.0479   0.0434 0.0576
## 3 2013-01-01 12:03:00 ARKANSAS_…  2.13 0.447  5.64 0.125  0.0499   0.0455 0.0625
## 4 2013-01-01 18:00:00 ARKANSAS_…  1.97 0.381  5.71 0.111  0.0426   0.0453 0.0555
## 5 2013-01-01 23:59:00 ARKANSAS_…  1.88 0.352  5.82 0.103  0.0393   0.0448 0.0516
## 6 2013-01-02 06:03:00 ARKANSAS_…  1.80 0.290  5.63 0.0959 0.0367   0.0443 0.0982
## # … with 7 more variables: NO3O2 <dbl>, PHYTO <dbl>, PON <dbl>, SOD_T <dbl>,
## #   TN <dbl>, TP <dbl>, WTEMP <dbl>


unique(Grab[c(4,6)])
## # A tibble: 52 × 2
##    Pcode    Units     
##    <chr>    <chr>     
##  1 COND     uS/cm @25C
##  2 DOSAT    %         
##  3 DO       mg/l      
##  4 PH       std units 
##  5 NO3_N    mg/l      
##  6 TEMPWC   degC      
##  7 FLOW_CMS cms       
##  8 TSS      mg/l      
##  9 TURB     NTU       
## 10 NO3O2_N  mg/l      
## # … with 42 more rows


library(readr)
Base_NoFrank <- read_delim("Base_Epa_Model_No_Franklin_Stp_Discharge.txt", 
    delim = "\t", escape_double = FALSE, 
    col_types = cols(`Date-Time` = col_datetime(format = "%m/%d/%Y %H:%M"), 
        CBOD = col_number(), DIP = col_number(), 
        DO = col_number(), DON = col_number(), 
        DOP = col_number(), FLOW_CMS = col_number(), 
        NH3_N = col_number(), NO3O2_N = col_number(), 
        PHYTO = col_number(), PON = col_number(), 
        SOD_T = col_number(), TN = col_number(), 
        TP = col_number(), WTEMP = col_number()), 
    trim_ws = TRUE)
head(Base_NoFrank)
## # A tibble: 6 × 16
##   `Date-Time`         Station_ID  CBOD   DIP    DO    DON    DOP FLOW_CMS  NH3_N
##   <dttm>              <chr>      <dbl> <dbl> <dbl>  <dbl>  <dbl>    <dbl>  <dbl>
## 1 2011-01-01 00:00:00 ARKANSAS_… 0     0      0    0      0         0.210 0     
## 2 2011-01-02 00:04:00 ARKANSAS_… 0.742 0.104  5.42 0.0231 0.0115    0.204 0.0116
## 3 2011-01-02 23:58:00 ARKANSAS_… 0.742 0.104  5.35 0.0231 0.0115    0.204 0.0646
## 4 2011-01-04 00:01:00 ARKANSAS_… 0.742 0.104 10.7  0.0231 0.0115    0.204 0.0633
## 5 2011-01-05 00:05:00 ARKANSAS_… 0.742 0.104 10.7  0.0231 0.0115    0.204 0.0621
## 6 2011-01-05 23:58:00 ARKANSAS_… 0.742 0.104 10.7  0.0231 0.0115    0.204 0.0610
## # … with 7 more variables: NO3O2_N <dbl>, PHYTO <dbl>, PON <dbl>, SOD_T <dbl>,
## #   TN <dbl>, TP <dbl>, WTEMP <dbl>



library(readr)
NPDES <- read_delim("Point_Source_Wastestream_Dataset.txt", 
    delim = "\t", escape_double = FALSE, 
    col_types = cols(`Date-Time` = col_datetime(format = "%m/%d/%Y %H:%M"), 
        CBOD = col_number(), DO = col_number(), 
        FLOW = col_number(), NH3 = col_number(), 
        NOX = col_number(), `NP RATIO` = col_number(), 
        ORGN = col_number(), ORGP = col_number(), 
        PO4 = col_number(), TKN = col_number(), 
        TN = col_number(), TP = col_number(), 
        TSS = col_number(), WTEMP = col_number()), 
    trim_ws = TRUE)

Grab data table contains grab sample results collected from 171 locations throughout the Harpeth River watershed


#Subset Data into Logical WASP Segments
# HARPETH_1 = Mouth of the Harpeth
# JONES_167 = Jones Creek
# HARPETH_2 = Main Stem w/o Jones Creek
# TURNBULL_143 = Turnbull Creek
# S_HARPETH_124 = South Harpeth
# HARPETH_34 = West Harpeth
# HARPETH_62 = Upstream of Franklin STP
# LSPC188PERO = Downstream of Franklin STP
# LSPC183PERO = Headwaters
# LSPC192PERO = Hwy100
# LSPC197PERO = 70s Bridge

sites=c("HARPETH_1","JONES_167","HARPETH_2", "TURNBULL_143", "S_HARPETH_124", "HARPETH_34", "LSPC188PERO", "HARPETH_62", "LSPC183PERO","LSPC192PERO", "LSPC197PERO")

colnames(BaseModel)[1] <- "Date"

BaseModel$TPLoad<-BaseModel$TP*0.00220462*BaseModel$FLOW_CMS*86400

BaseModel$TPLoad.csum <- ave(BaseModel$TPLoad, BaseModel$Station_ID, FUN=cumsum)

BaseModel$NP<-BaseModel$TN/BaseModel$TP

BaseModel1<-BaseModel

BaseModel<-BaseModel[BaseModel$Station_ID==sites,]

colnames(Grab)[3] <- "Date"

colnames(Cont)[2] <- "Date"

colnames(Natural)[1] <- "Date"

colnames(Base_PS)[1] <- "Date"

colnames(NPDES)[1] <- "Date"

Natural$NP<-Natural$TN/Natural$TP

Natural$TPLoad<-Natural$TP*0.00220462*Natural$FLOW_CMS*86400

Natural$TPLoad.csum <- ave(Natural$TPLoad, Natural$Station_ID, FUN=cumsum)

Natural1<-Natural

Natural<-Natural[Natural$Station_ID==sites,]



Base_PS<-Base_PS[Base_PS$Station_ID==sites,]

Base_PS$TPLoad<-Base_PS$TP*0.00220462*Base_PS$FLOW_CMS*86400

Base_PS1<-Base_PS

colnames(Base_NoFrank)[1] <- "Date"

Base_NoFrank$TPLoad<-Base_NoFrank$TP*0.00220462*Base_NoFrank$FLOW_CMS*86400

Base_NoFrank$TPLoad.csum <- ave(Base_NoFrank$TPLoad, Base_NoFrank$Station_ID, FUN=cumsum)

Base_NoFrank$NP<-Base_NoFrank$TN/Base_NoFrank$TP

Base_NoFrank1<-Base_NoFrank

Base_NoFrank<-Base_NoFrank[Base_NoFrank$Station_ID==sites,]


NPDES$TPLoad<-NPDES$TP*0.00220462*NPDES$FLOW*86400

NPDES$TPLoad.csum <- ave(NPDES$TPLoad, NPDES$Station_ID, FUN=cumsum)

NPDES$NP<-NPDES$TN/NPDES$TP

NPDES1<-NPDES




#####################################################################################################

molten.data <- melt(Grab, id = c("Date","Station_ID","Pcode","Units","DataSource"))

Grab.w<-dcast(molten.data, Date+Station_ID~Pcode, mean)

Grab.w$NP<-Grab.w$TN/Grab.w$TP

head(Grab.w, 10)
##                   Date             Station_ID ALK BOD30 BOD5 CBOD30 CBOD5
## 1  2012-01-03 08:33:00 11NPSWRD_WQX_NATR_BUBR NaN   NaN  NaN    NaN   NaN
## 2  2012-01-03 09:35:00           TNW000003300 NaN   NaN  NaN    NaN   NaN
## 3  2012-01-03 10:45:00           TNW000002790 NaN   NaN  NaN    NaN     1
## 4  2012-01-03 11:15:00           TNW000003508 NaN   NaN  NaN    NaN   NaN
## 5  2012-01-03 11:55:00           TNW000002791 NaN   NaN  NaN    NaN   NaN
## 6  2012-01-03 12:20:00           TNW000005654 NaN   NaN  NaN    NaN   NaN
## 7  2012-01-03 13:15:00           TNW000006704 NaN   NaN  NaN    NaN   NaN
## 8  2012-01-04 08:30:00           TNW000003576 NaN   NaN  NaN    NaN     1
## 9  2012-01-04 09:43:00           TNW000006443 NaN   NaN  NaN    NaN   NaN
## 10 2012-01-04 10:25:00           TNW000003306 NaN   NaN  NaN    NaN   NaN
##    CBOD90 CHLA COND    DO DOSAT FLOW_CFS   FLOW_CMS NBOD90  NH3_N NO2 NO2_N NO3
## 1     NaN  NaN  306 12.55  90.3      NaN 0.02123764    NaN    NaN NaN   NaN NaN
## 2     NaN  NaN  295 14.36   NaN      NaN        NaN    NaN    NaN NaN   NaN NaN
## 3     NaN  NaN  373 12.81   NaN      NaN        NaN    NaN 0.0165 NaN   NaN NaN
## 4     NaN  NaN  242 15.58   NaN      NaN        NaN    NaN    NaN NaN   NaN NaN
## 5     NaN  NaN  376 13.39   NaN      NaN        NaN    NaN    NaN NaN   NaN NaN
## 6     NaN  NaN  246 15.08   NaN      NaN        NaN    NaN    NaN NaN   NaN NaN
## 7     NaN  NaN  234 16.51   NaN      NaN        NaN    NaN    NaN NaN   NaN NaN
## 8     NaN  NaN  544 13.33   NaN      NaN        NaN    NaN 0.0165 NaN   NaN NaN
## 9     NaN  NaN  317 14.70   NaN      NaN        NaN    NaN 0.0165 NaN   NaN NaN
## 10    NaN  NaN  421 17.50   NaN      NaN        NaN    NaN 0.0165 NaN   NaN NaN
##    NO3_N NO3O2_N OPO4 OPO4_P ORGN ORP PERI%N PERI%P PERI_AI PERI_ASH PERI_NP
## 1   3.03     NaN  NaN    NaN  NaN NaN    NaN    NaN     NaN      NaN     NaN
## 2    NaN    0.43  NaN    NaN  NaN NaN    NaN    NaN     NaN      NaN     NaN
## 3    NaN    1.10  NaN    NaN  NaN NaN    NaN    NaN     NaN      NaN     NaN
## 4    NaN    0.15  NaN    NaN  NaN NaN    NaN    NaN     NaN      NaN     NaN
## 5    NaN    1.10  NaN    NaN  NaN NaN    NaN    NaN     NaN      NaN     NaN
## 6    NaN    0.41  NaN    NaN  NaN NaN    NaN    NaN     NaN      NaN     NaN
## 7    NaN    0.76  NaN    NaN  NaN NaN    NaN    NaN     NaN      NaN     NaN
## 8    NaN    1.30  NaN    NaN  NaN NaN    NaN    NaN     NaN      NaN     NaN
## 9    NaN    1.70  NaN    NaN  NaN NaN    NaN    NaN     NaN      NaN     NaN
## 10   NaN    0.92  NaN    NaN  NaN NaN    NaN    NaN     NaN      NaN     NaN
##    PERI_TN PERI_TP PERIAFDM PERICHLA   PH SOD_20 SOD_T SRP TBOD90 TDP TDS
## 1      NaN     NaN      NaN      NaN 7.76    NaN   NaN NaN    NaN NaN NaN
## 2      NaN     NaN      NaN      NaN 7.85    NaN   NaN NaN    NaN NaN NaN
## 3      NaN     NaN      NaN      NaN 8.04    NaN   NaN NaN    NaN NaN NaN
## 4      NaN     NaN      NaN      NaN 8.27    NaN   NaN NaN    NaN NaN NaN
## 5      NaN     NaN      NaN      NaN 8.13    NaN   NaN NaN    NaN NaN NaN
## 6      NaN     NaN      NaN      NaN 8.27    NaN   NaN NaN    NaN NaN NaN
## 7      NaN     NaN      NaN      NaN 8.45    NaN   NaN NaN    NaN NaN NaN
## 8      NaN     NaN      NaN      NaN 7.41    NaN   NaN NaN    NaN NaN NaN
## 9      NaN     NaN      NaN      NaN 8.31    NaN   NaN NaN    NaN NaN NaN
## 10     NaN     NaN      NaN      NaN 8.59    NaN   NaN NaN    NaN NaN NaN
##    TEMPWC   TKN  TN    TP  TS TSS TURB  NP
## 1    0.40   NaN NaN   NaN NaN NaN  NaN NaN
## 2    4.33 0.065 NaN 0.006 NaN   5 0.68 NaN
## 3    6.16 0.065 NaN 0.210 NaN   5 3.07 NaN
## 4    2.61 0.065 NaN 0.013 NaN   5 0.32 NaN
## 5    6.36 0.200 NaN 0.210 NaN   5 2.81 NaN
## 6    6.19 0.065 NaN 0.044 NaN   5 0.74 NaN
## 7    3.39 0.065 NaN 0.082 NaN   5 1.75 NaN
## 8    4.62 0.065 NaN 0.230 NaN   5 1.36 NaN
## 9    4.60 0.065 NaN 0.210 NaN   5 0.52 NaN
## 10   5.46 0.200 NaN 0.033 NaN   5 0.88 NaN



molten.cont <- melt(Cont, id = c("Date","Station_ID","Pcode","Units"))

Cont.w<-dcast(molten.cont, Date+Station_ID~Pcode, mean)

head(Cont.w, 10)
##          Date    Station_ID COND  DO DOSAT FLOW_CFS    FLOW_CMS  PH TEMPWC TURB
## 1  2011-01-01 USGS_03432350  NaN NaN   NaN   828.00 23.44634932 NaN    NaN  NaN
## 2  2011-01-01 USGS_03432400  NaN NaN   NaN   919.00 26.02318239 NaN    NaN  NaN
## 3  2011-01-01 USGS_03433500  NaN NaN   NaN   456.00 12.91248223 NaN    NaN  NaN
## 4  2011-01-01 USGS_03433640  NaN NaN   NaN     3.83  0.10845352 NaN    NaN  NaN
## 5  2011-01-01 USGS_03434500  NaN NaN   NaN   464.00 13.13901701 NaN    NaN  NaN
## 6  2011-01-02 USGS_03432350  NaN NaN   NaN  1020.00 28.88318394 NaN    NaN  NaN
## 7  2011-01-02 USGS_03432400  NaN NaN   NaN  1330.00 37.66140651 NaN    NaN  NaN
## 8  2011-01-02 USGS_03433500  NaN NaN   NaN  2220.00 62.86340034 NaN    NaN  NaN
## 9  2011-01-02 USGS_03433640  NaN NaN   NaN     1.44  0.04077626 NaN    NaN  NaN
## 10 2011-01-02 USGS_03434500  NaN NaN   NaN  2130.00 60.31488411 NaN    NaN  NaN

#Data Summary

Subset of 11 modeled segments in WASP


ggplot(data = BaseModel, aes(x = Date, y = FLOW_CMS, group = Station_ID, colour = Station_ID)) +
    geom_line() +
    facet_wrap(~ Station_ID)







ggplot(data = BaseModel, aes(x = Date, y = TPLoad, group = Station_ID, colour = Station_ID, )) +
  ylab("TP Load (lbs/day)") +
    geom_line() +
    facet_wrap(~ Station_ID)



#Boxplot - look at data distribution



tp<-ggplot(BaseModel, aes(x = Station_ID, y = TP, fill = Station_ID)) + 
  geom_boxplot(outlier.colour="black", outlier.shape=16, outlier.size=2) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+
  ylab("TP (mg/L)")

tp



tn<-ggplot(BaseModel, aes(x = Station_ID, y = TN, fill = Station_ID)) + 
  geom_boxplot(outlier.colour="black", outlier.shape=16, outlier.size=2) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+
  ylab("TN (mg/L)")

tn



np<-ggplot(BaseModel, aes(x = Station_ID, y = NP, fill = Station_ID)) + 
  geom_boxplot(outlier.colour="black", outlier.shape=16, outlier.size=2) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+
  ylab("TN:TP Ratio")

np


LSPC188_tp<-ggplot(BaseModel,aes(x=Date,y=TP, fill="Base Model"))+geom_line(data=subset(BaseModel,Station_ID=="LSPC188PERO")) +theme_bw()+
  labs(y="TP (mg/L)")+ylim(min(BaseModel$TP),max(BaseModel$TP))+
  theme(
    axis.text.y=element_text(size=11),axis.title.y=element_text(size=12, vjust=1.2),
    axis.text.x=element_blank(),axis.title.x=element_blank(), 
    panel.border = element_rect(size=.8, colour = "black"))+
  geom_point(data=subset(Grab.w, Station_ID==c("Frank_Down","Frank_CottonLn","Franklin_Site_2","Franklin_Site_3")), aes(x=Date, y=TP, color=Station_ID))+
  geom_line(data=subset(Natural, Station_ID=="LSPC188PERO"), aes(x=Date, y=TP, color="Natural"))

LSPC188_tp


LSPC188_tn<-ggplot(BaseModel,aes(x=Date,y=TN))+geom_line(data=subset(BaseModel,Station_ID=="LSPC188PERO")) +theme_bw()+
  labs(y="TN (mg/L)")+ylim(min(BaseModel$TN),max(BaseModel$TN))+
  theme(
    axis.text.y=element_text(size=11),axis.title.y=element_text(size=12, vjust=1.2),
    axis.text.x=element_blank(),axis.title.x=element_blank(), 
    panel.border = element_rect(size=.8, colour = "black"))+
  geom_point(data=subset(Grab.w, Station_ID==c("Frank_Down","Frank_CottonLn","Franklin_Site_2","Franklin_Site_3")), aes(x=Date, y=TN, color=Station_ID))+
  geom_line(data=subset(Natural, Station_ID=="LSPC188PERO"), aes(x=Date, y=TN, color="Natural"))

LSPC188_tn


LSPC188_flow<-ggplot(BaseModel,aes(x=Date,y=FLOW_CMS))+geom_line(data=subset(BaseModel,Station_ID=="LSPC188PERO")) +theme_bw()+
  labs(y="Flow (cms)")+ylim(min(BaseModel$FLOW_CMS),200)+
  theme(
    axis.text.y=element_text(size=11),axis.title.y=element_text(size=12, vjust=1.2),
    axis.text.x=element_blank(),axis.title.x=element_blank(), 
    panel.border = element_rect(size=.8, colour = "black"))+
  geom_point(data=subset(Cont.w, Station_ID==c("USGS_03432400")), aes(x=Date, y=FLOW_CMS, color=Station_ID))+
  geom_line(data=subset(Natural, Station_ID=="LSPC188PERO"), aes(x=Date, y=FLOW_CMS, color="Natural"))

LSPC188_flow




LSPC188_np<-ggplot(BaseModel,aes(x=Date,y=NP))+geom_line(data=subset(BaseModel,Station_ID=="LSPC188PERO")) +theme_bw()+
  labs(y="TN:TP ratio", x="Date")+ylim(min(BaseModel$NP),max(BaseModel$NP))+
  theme(
    axis.text.y=element_text(size=11),axis.title.y=element_text(size=12, vjust=1.2),
    axis.text.x=element_text(size=11),axis.title.x=element_text(size=11), 
    panel.border = element_rect(size=.8, colour = "black"))+
  geom_point(data=subset(Grab.w, Station_ID==c("Frank_Down","Frank_CottonLn","Franklin_Site_2","Franklin_Site_3")), aes(x=Date, y=NP, color=Station_ID))+
  geom_line(data=subset(Natural, Station_ID=="LSPC188PERO"), aes(x=Date, y=NP, color="Natural"))

LSPC188_np


grid.arrange(LSPC188_tp,LSPC188_tn,LSPC188_flow,LSPC188_np, ncol=1,padding=100)



HARPETH_62<-ggplot(BaseModel,aes(x=Date,y=TP))+geom_line(data=subset(BaseModel,Station_ID=="HARPETH_62")) +theme_bw()+
  labs(y="TP (mg/L)")+ylim(min(BaseModel$TP),max(BaseModel$TP))+
  theme(
    axis.text.y=element_text(size=16),axis.title.y=element_text(size=12, vjust=1.2),
    axis.text.x=element_blank(),axis.title.x=element_blank(), 
    panel.border = element_rect(size=.8, colour = "black"))+
  geom_point(data=subset(Grab.w, Station_ID==c("Frank_Up","Frank_Eff","Franklin_Site_1")), aes(x=Date, y=TP, color=Station_ID))


#HARPETH_62







#ppi=300
#png("awesome stacked time series.png",width=12*ppi, height=8*ppi, res=ppi)
#grid.arrange(LSPC188_tp,HARPETH_62, ncol=1)
#dev.off()

LSPC188_DO<-ggplot(BaseModel,aes(x=Date,y=DO, fill="Base Model"))+geom_line(data=subset(BaseModel,Station_ID=="LSPC188PERO")) +theme_bw()+
  labs(y="DO (mg/L)")+ylim(min(BaseModel$DO),max(BaseModel$DO))+
  theme(
    axis.text.y=element_text(size=11),axis.title.y=element_text(size=12, vjust=1.2),
    axis.text.x=element_text(size=11),axis.title.x=element_text(size=11), 
    panel.border = element_rect(size=.8, colour = "black"))+
  geom_point(data=subset(Grab.w, Station_ID==c("Frank_Down","Frank_CottonLn","Franklin_Site_2","Franklin_Site_3")), aes(x=Date, y=DO, color=Station_ID))+
  geom_line(data=subset(Natural, Station_ID=="LSPC188PERO"), aes(x=Date, y=DO, colour="Natural"))

LSPC188_DO


LSPC188_DO_PSinc<-ggplot(Base_PS,aes(x=Date,y=DO,fill="Base Model x 2.25 PS"))+geom_line(data=subset(Base_PS,Station_ID=="LSPC188PERO")) +theme_bw()+
  labs(y="DO (mg/L)")+ylim(min(Base_PS$DO),max(Base_PS$DO))+
  theme(
    axis.text.y=element_text(size=11),axis.title.y=element_text(size=12, vjust=1.2),
    axis.text.x=element_text(size=11),axis.title.x=element_text(size=11), 
    panel.border = element_rect(size=.8, colour = "black"))+
  geom_point(data=subset(Grab.w, Station_ID==c("Frank_Down","Frank_CottonLn","Franklin_Site_2","Franklin_Site_3")), aes(x=Date, y=DO, colour=Station_ID))+
  geom_line(data=subset(Natural, Station_ID=="LSPC188PERO"), aes(x=Date, y=DO, color="Natural"))

LSPC188_DO_PSinc


LSPC188_flow<-ggplot(BaseModel,aes(x=Date,y=FLOW_CMS))+geom_line(data=subset(BaseModel,Station_ID=="LSPC188PERO")) +theme_bw()+
  labs(y="Flow (cms)")+ylim(min(BaseModel$FLOW_CMS),200)+
  theme(
    axis.text.y=element_text(size=11),axis.title.y=element_text(size=12, vjust=1.2),
    axis.text.x=element_text(size=11),axis.title.x=element_text(size=11), 
    panel.border = element_rect(size=.8, colour = "black"))+
  geom_point(data=subset(Cont.w, Station_ID==c("USGS_03432400")), aes(x=Date, y=FLOW_CMS, color=Station_ID))+
  geom_line(data=subset(Natural, Station_ID=="LSPC188PERO"), aes(x=Date, y=FLOW_CMS, color="Natural"))

LSPC188_flow



grid.arrange(LSPC188_DO,LSPC188_DO_PSinc,LSPC188_flow, ncol=1,padding=100)




p<-ggplot(BaseModel, aes(x=Date, y=DO, fill=Station_ID)) + ggtitle("Base Model/Current Conditions") + theme_bw() +
  geom_line(aes(colour=Station_ID, group=Station_ID)) #+
  #geom_point(data=subset(Cont.w, Station_ID==c("USGS_03434500","HARPE062")) , aes(x=Date, y=DO))

p



p1<-ggplot(Natural, aes(x=Date, y=DO, fill=Station_ID)) + ggtitle("Natural Conditions") + theme_bw() +
  geom_line(aes(colour=Station_ID, group=Station_ID))

p1



p3<-ggplot(Base_PS, aes(x=Date, y=DO, fill=Station_ID)) + theme_bw() +
  geom_line(aes(colour=Station_ID, group=Station_ID))

p3



monthOrder <- c('Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec')
BaseModel$Month <- factor(format(BaseModel$Date, "%b"), levels = monthOrder)
BaseModel$Year <- factor(format(BaseModel$Date, "%Y"))
ggplot(BaseModel, aes(Month, TPLoad)) + geom_boxplot(data=subset(BaseModel,Station_ID=="LSPC188PERO")) + stat_boxplot(geom ='errorbar') + ggtitle("TP Load (lbs/day)")


ggplot(BaseModel,aes(Month,TPLoad)) + geom_bar(stat="identity") + ggtitle("TP Load (lbs/month)")



Natural$Month <- factor(format(Natural$Date, "%b"), levels = monthOrder)
Natural$Year <- factor(format(Natural$Date, "%Y"))
ggplot(Natural, aes(Month, TPLoad)) + geom_boxplot(data=subset(Natural,Station_ID=="LSPC188PERO")) + stat_boxplot(geom ='errorbar') + ggtitle("Natural Model - TP Load (lbs/day)")



Base_NoFrank1$Month <- factor(format(Base_NoFrank1$Date, "%b"), levels = monthOrder)
Base_NoFrank1$Year <- factor(format(Base_NoFrank1$Date, "%Y"))
ggplot(Base_NoFrank1, aes(Month, TPLoad)) + geom_boxplot(data=subset(Base_NoFrank1,Station_ID=="LSPC188PERO")) + stat_boxplot(geom ='errorbar') + ggtitle("Base_NoFrank1 Model - TP Load (lbs/day)")



Base_PS1$Month <- factor(format(Base_PS1$Date, "%b"), levels = monthOrder)
Base_PS1$Year <- factor(format(Base_PS1$Date, "%Y"))
ggplot(Base_PS1, aes(Month, TPLoad)) + geom_boxplot(data=subset(Base_PS1,Station_ID=="LSPC188PERO")) + stat_boxplot(geom ='errorbar') + ggtitle("Base_PS1 Model - TP Load (lbs/day)")


##################

#LSPC188_tpload <- subset(BaseModel, Station_ID==c("LSPC188PERO","HARPETH_62"), select=c(Station_ID, Year, Month, TPLoad))
#meanTPLoad_62 <- subset(BaseModel, Station_ID=="HARPETH_62", select=c(Station_ID, Year, Month, TPLoad))
#meanTPLoad_188 <- subset(BaseModel, Station_ID=="LSPC188PERO", select=c(Station_ID, Year, Month, TPLoad))

#library(plyr)
#tpload62<-ddply(meanTPLoad_62, c("Month","Station_ID"), summarise, x = mean(TPLoad))

#meanTPLoad_62 <- join(meanTPLoad_62, tpload62, match="all")






#ggplot(LSPC188_tpload,aes(Year,TPLoad,colour=Station_ID)) +
#geom_point(data=LSPC188_tpload,size=I(2),alpha=I(0.6)) + 
#geom_line(data=meanTPLoad_62, aes(Month,x), size=I(1.5),alpha=I(0.6)) + 
##geom_line(data=mean(meanTPLoad_62$TPLoad),size=I(1.5),alpha=I(0.4)) + 
#theme_grey(base_size=15) +
#theme(legend.title = element_blank(), legend.position=c(.85,.85), axis.title.y=element_blank(),axis.text.x=element_blank()) + 
#ggtitle("TP Load (lbs/day)") + facet_grid(. ~ Month) + 
#xlab(paste("Years: 2011 to 2019"))













#theme_set(theme_classic())

#BaseModel$Date.ts<-as.ts(BaseModel$Date)

# Subset data
#Base_LSPC188 <- window(subset(BaseModel,Station_ID=="LSPC188PERO"), start=c(2012, 1), end=c(2019, 8))  # subset a smaller timewindow

# Plot
#ggseasonplot(subset(BaseModel,Station_ID=="LSPC188PERO"), x=BaseModel$Date.ts) + labs(title="Seasonal plot: International Airline Passengers")
#ggseasonplot(subset(BaseModel,Station_ID=="LSPC188PERO")) + labs(title="Seasonal plot: Air temperatures at Nottingham Castle")
p <- ggplot(data=BaseModel1, aes(x=Date, y=DO, color=Station_ID)) +
  geom_line(data=subset(BaseModel1,Station_ID=="LSPC188PERO")) + 
  xlab("Year") +
  geom_line(data=subset(BaseModel1, Station_ID==sites))

p

p <- ggplot(data=Base_PS1, aes(x=Date, y=DO, color=Station_ID)) +
  geom_line(data=subset(Base_PS1,Station_ID=="LSPC188PERO")) + 
  xlab("Year") +
  geom_line(data=subset(Base_PS1, Station_ID==sites))

p

p <- ggplot(data=Natural1, aes(x=Date, y=DO, color=Station_ID)) +
  geom_line(data=subset(Natural1,Station_ID=="LSPC188PERO")) + 
  xlab("Year") +
  geom_line(data=subset(Natural1, Station_ID==sites))

p

LSPC188.base<-subset(BaseModel1,Station_ID=="LSPC188PERO")
dy.tpload <- xts(x = LSPC188.base$TPLoad, order.by = LSPC188.base$Date)

LSPC188.natural<-subset(Natural1,Station_ID=="LSPC188PERO")
dy.Nattpload <- xts(x = LSPC188.natural$TPLoad, order.by = LSPC188.natural$Date)

LSPC188.tpload<-cbind(dy.tpload,dy.Nattpload)


LSPC188.tpload <- na.locf(LSPC188.tpload)

p <- dygraph(LSPC188.tpload, main="TP Load at LSPC 188 - Below Frank", ylab="TP Load (lbs/day)") %>%
  dySeries("dy.tpload", drawPoints = TRUE, strokePattern="dashed", color="red") %>%
  dySeries("dy.Nattpload", drawPoints= T, color="blue") %>%
  dyOptions(labelsUTC = TRUE, fillGraph=TRUE, fillAlpha=0.1, drawGrid = FALSE, colors = RColorBrewer::brewer.pal(3, "Set2")) %>%
  dyRangeSelector() %>%
  dyCrosshair(direction = "vertical") %>%
  dyHighlight(highlightCircleSize = 5, highlightSeriesBackgroundAlpha = 0.2, hideOnMouseOut = FALSE)  %>%
  dyRoller(rollPeriod = 1)


p
LSPC192.base<-subset(BaseModel1,Station_ID=="LSPC192PERO")
dy.DO <- xts(x = LSPC192.base$DO, order.by = LSPC192.base$Date)

LSPC192.natural<-subset(Natural1,Station_ID=="LSPC192PERO")
dy.NatDO <- xts(x = LSPC192.natural$DO, order.by = LSPC192.natural$Date)

LSPC192.Base_PS<-subset(Base_PS1,Station_ID=="LSPC192PERO")
dy.Base_PSDO <- xts(x = LSPC192.Base_PS$DO, order.by = LSPC192.Base_PS$Date)

LSPC192.Cont<-subset(Cont.w,Station_ID==c("USGS_03433500","HARPE062"))
dy.ContDO <- xts(x = LSPC192.Cont$DO, order.by = LSPC192.Cont$Date)

LSPC192.DO<-cbind(dy.DO,dy.NatDO,dy.Base_PSDO,dy.ContDO)


#LSPC192.DO <- na.locf(LSPC192.DO)


p1 <- dygraph(LSPC192.DO, main="Dissolved Oxygen at LSPC 192 - Hwy 100", ylab="DO (mg/L)") %>%
  dySeries("dy.DO", drawPoints = TRUE, strokePattern="dashed", color="red") %>%
  dySeries("dy.NatDO", drawPoints= T, color="blue") %>%
  dyOptions(labelsUTC = TRUE, fillGraph=TRUE, fillAlpha=0.1, drawGrid = FALSE, colors = RColorBrewer::brewer.pal(4, "Set1")) %>%
  dyRangeSelector() %>%
  dyCrosshair(direction = "vertical") %>%
  dyHighlight(highlightCircleSize = 5, highlightSeriesBackgroundAlpha = 0.2, hideOnMouseOut = FALSE)  %>%
  dyRoller(rollPeriod = 1)


p1
LSPC188.base<-subset(BaseModel1,Station_ID=="LSPC188PERO")
dy.DO <- xts(x = LSPC188.base$DO, order.by = LSPC188.base$Date)

LSPC188.natural<-subset(Natural1,Station_ID=="LSPC188PERO")
dy.NatDO <- xts(x = LSPC188.natural$DO, order.by = LSPC188.natural$Date)

LSPC188.Base_PS<-subset(Base_PS1,Station_ID=="LSPC188PERO")
dy.Base_PSDO <- xts(x = LSPC188.Base_PS$DO, order.by = LSPC188.Base_PS$Date)

LSPC188.Cont<-subset(Cont.w,Station_ID==c("USGS_03432400","USGS_03432350"))
dy.ContDO <- xts(x = LSPC188.Cont$DO, order.by = LSPC188.Cont$Date)

LSPC188.Grab<-subset(Grab.w,Station_ID==c("Frank_Down","Frank_CottonLn", "Frank_Eff", "Franklin_Site_2", "Franklin_Site_3"))
## Warning in Station_ID == c("Frank_Down", "Frank_CottonLn", "Frank_Eff", : longer
## object length is not a multiple of shorter object length
dy.GrabDO <- xts(x = LSPC188.Grab$DO, order.by = LSPC188.Grab$Date)


LSPC188.DO<-cbind(dy.DO,dy.NatDO,dy.Base_PSDO,dy.ContDO,dy.GrabDO)


#LSPC188.DO <- na.locf(LSPC188.DO)






p2 <- dygraph(LSPC188.DO, main="Dissolved Oxygen at LSPC 188 - Below Franklin STP", ylab="DO (mg/L)") %>%
  dySeries("dy.DO", drawPoints = TRUE, strokePattern="dashed", color="red") %>%
  dySeries("dy.NatDO", drawPoints= T, color="blue") %>%
  dyOptions(labelsUTC = TRUE, fillGraph=TRUE, fillAlpha=0.1, drawGrid = FALSE, colors = RColorBrewer::brewer.pal(5, "Set1")) %>%
  dyRangeSelector() %>%
  dyCrosshair(direction = "vertical") %>%
  dyHighlight(highlightCircleSize = 5, highlightSeriesBackgroundAlpha = 0.2, hideOnMouseOut = FALSE)  %>%
  dyRoller(rollPeriod = 1)


p2
HARPETH_62.base<-subset(BaseModel1,Station_ID=="HARPETH_62")
dy.DO <- xts(x = HARPETH_62.base$DO, order.by = HARPETH_62.base$Date)

HARPETH_62.natural<-subset(Natural1,Station_ID=="HARPETH_62")
dy.NatDO <- xts(x = HARPETH_62.natural$DO, order.by = HARPETH_62.natural$Date)

HARPETH_62.Base_PS<-subset(Base_PS1,Station_ID=="HARPETH_62")
dy.Base_PSDO <- xts(x = HARPETH_62.Base_PS$DO, order.by = HARPETH_62.Base_PS$Date)

HARPETH_62.Cont<-subset(Cont.w,Station_ID==c("usgs_03432350","USGS_0343233905"))
dy.ContDO <- xts(x = HARPETH_62.Cont$DO, order.by = HARPETH_62.Cont$Date)

HARPETH_62.Grab<-subset(Grab.w,Station_ID==c("Frank_Up","Franklin_Site_1"))
## Warning in Station_ID == c("Frank_Up", "Franklin_Site_1"): longer object length
## is not a multiple of shorter object length
dy.GrabDO <- xts(x = HARPETH_62.Grab$DO, order.by = HARPETH_62.Grab$Date)


HARPETH_62.DO<-cbind(dy.DO,dy.NatDO,dy.Base_PSDO,dy.ContDO,dy.GrabDO)


#HARPETH_62.DO <- na.locf(HARPETH_62.DO)






p3 <- dygraph(HARPETH_62.DO, main="Dissolved Oxygen at Harpeth 62 - Upstream of Franklin STP", ylab="DO (mg/L)") %>%
  dySeries("dy.DO", drawPoints = TRUE, strokePattern="dashed", color="red") %>%
  dySeries("dy.NatDO", drawPoints= T, color="blue") %>%
  dyOptions(labelsUTC = TRUE, fillGraph=TRUE, fillAlpha=0.1, drawGrid = FALSE, colors = RColorBrewer::brewer.pal(5, "Set1")) %>%
  dyRangeSelector() %>%
  dyCrosshair(direction = "vertical") %>%
  dyHighlight(highlightCircleSize = 5, highlightSeriesBackgroundAlpha = 0.2, hideOnMouseOut = FALSE)  %>%
  dyRoller(rollPeriod = 1)


p3
JONES_167.base<-subset(BaseModel1,Station_ID=="JONES_167")
dy.DO <- xts(x = JONES_167.base$DO, order.by = JONES_167.base$Date)

JONES_167.natural<-subset(Natural1,Station_ID=="JONES_167")
dy.NatDO <- xts(x = JONES_167.natural$DO, order.by = JONES_167.natural$Date)

JONES_167.Base_PS<-subset(Base_PS1,Station_ID=="JONES_167")
dy.Base_PSDO <- xts(x = JONES_167.Base_PS$DO, order.by = JONES_167.Base_PS$Date)

#JONES_167.Cont<-subset(Cont.w,Station_ID==c("usgs_03432350","USGS_0343233905"))
#dy.ContDO <- xts(x = JONES_167.Cont$DO, order.by = JONES_167.Cont$Date)

#JONES_167.Grab<-subset(Grab.w,Station_ID==c("Frank_Up","Franklin_Site_1"))
#dy.GrabDO <- xts(x = JONES_167.Grab$DO, order.by = JONES_167.Grab$Date)


JONES_167.DO<-cbind(dy.DO,dy.NatDO,dy.Base_PSDO)


#JONES_167.DO <- na.locf(JONES_167.DO)






p4 <- dygraph(JONES_167.DO, main="Dissolved Oxygen at JONES_167 - Mouth of Jones Creek", ylab="DO (mg/L)") %>%
  dySeries("dy.DO", drawPoints = TRUE, strokePattern="dashed", color="red") %>%
  dySeries("dy.NatDO", drawPoints= T, color="blue") %>%
  dyOptions(labelsUTC = TRUE, fillGraph=TRUE, fillAlpha=0.1, drawGrid = FALSE, colors = RColorBrewer::brewer.pal(5, "Set1")) %>%
  dyRangeSelector() %>%
  dyCrosshair(direction = "vertical") %>%
  dyHighlight(highlightCircleSize = 5, highlightSeriesBackgroundAlpha = 0.2, hideOnMouseOut = FALSE)  %>%
  dyRoller(rollPeriod = 1)


p4
LSPC183PERO.base<-subset(BaseModel1,Station_ID=="LSPC183PERO")
dy.DO <- xts(x = LSPC183PERO.base$DO, order.by = LSPC183PERO.base$Date)

LSPC183PERO.natural<-subset(Natural1,Station_ID=="LSPC183PERO")
dy.NatDO <- xts(x = LSPC183PERO.natural$DO, order.by = LSPC183PERO.natural$Date)

LSPC183PERO.Base_PS<-subset(Base_PS1,Station_ID=="LSPC183PERO")
dy.Base_PSDO <- xts(x = LSPC183PERO.Base_PS$DO, order.by = LSPC183PERO.Base_PS$Date)

#LSPC183PERO.Cont<-subset(Cont.w,Station_ID==c("usgs_03432350","USGS_0343233905"))
#dy.ContDO <- xts(x = LSPC183PERO.Cont$DO, order.by = LSPC183PERO.Cont$Date)

#LSPC183PERO.Grab<-subset(Grab.w,Station_ID==c("Frank_Up","Franklin_Site_1"))
#dy.GrabDO <- xts(x = LSPC183PERO.Grab$DO, order.by = LSPC183PERO.Grab$Date)


LSPC183PERO.DO<-cbind(dy.DO,dy.NatDO,dy.Base_PSDO)


#LSPC183PERO.DO <- na.locf(LSPC183PERO.DO)






p4 <- dygraph(LSPC183PERO.DO, main="Dissolved Oxygen at LSPC183PERO - Headwaters", ylab="DO (mg/L)") %>%
  dySeries("dy.DO", drawPoints = TRUE, strokePattern="dashed", color="red") %>%
  dySeries("dy.NatDO", drawPoints= T, color="blue") %>%
  dyOptions(labelsUTC = TRUE, fillGraph=TRUE, fillAlpha=0.1, drawGrid = FALSE, colors = RColorBrewer::brewer.pal(5, "Set1")) %>%
  dyRangeSelector() %>%
  dyCrosshair(direction = "vertical") %>%
  dyHighlight(highlightCircleSize = 5, highlightSeriesBackgroundAlpha = 0.2, hideOnMouseOut = FALSE)  %>%
  dyRoller(rollPeriod = 1)


p4
S_HARPETH_124.base<-subset(BaseModel1,Station_ID=="S_HARPETH_124")
dy.DO <- xts(x = S_HARPETH_124.base$DO, order.by = S_HARPETH_124.base$Date)

S_HARPETH_124.natural<-subset(Natural1,Station_ID=="S_HARPETH_124")
dy.NatDO <- xts(x = S_HARPETH_124.natural$DO, order.by = S_HARPETH_124.natural$Date)

S_HARPETH_124.Base_PS<-subset(Base_PS1,Station_ID=="S_HARPETH_124")
dy.Base_PSDO <- xts(x = S_HARPETH_124.Base_PS$DO, order.by = S_HARPETH_124.Base_PS$Date)

#S_HARPETH_124.Cont<-subset(Cont.w,Station_ID==c("usgs_03432350","USGS_0343233905"))
#dy.ContDO <- xts(x = S_HARPETH_124.Cont$DO, order.by = S_HARPETH_124.Cont$Date)

#S_HARPETH_124.Grab<-subset(Grab.w,Station_ID==c("Frank_Up","Franklin_Site_1"))
#dy.GrabDO <- xts(x = S_HARPETH_124.Grab$DO, order.by = S_HARPETH_124.Grab$Date)


S_HARPETH_124.DO<-cbind(dy.DO,dy.NatDO,dy.Base_PSDO)


#S_HARPETH_124.DO <- na.locf(S_HARPETH_124.DO)






p5 <- dygraph(S_HARPETH_124.DO, main="Dissolved Oxygen at S_HARPETH_124 - Mouth of South Harpeth", ylab="DO (mg/L)") %>%
  dySeries("dy.DO", drawPoints = TRUE, strokePattern="dashed", color="red") %>%
  dySeries("dy.NatDO", drawPoints= T, color="blue") %>%
  dyOptions(labelsUTC = TRUE, fillGraph=TRUE, fillAlpha=0.1, drawGrid = FALSE, colors = RColorBrewer::brewer.pal(5, "Set1")) %>%
  dyRangeSelector() %>%
  dyCrosshair(direction = "vertical") %>%
  dyHighlight(highlightCircleSize = 5, highlightSeriesBackgroundAlpha = 0.2, hideOnMouseOut = FALSE)  %>%
  dyRoller(rollPeriod = 1)


p5
DF <- melt(BaseModel, id=c("Date","Station_ID"))


ggplot(BaseModel, aes(x = TP, col = Station_ID))+
   stat_ecdf(lwd = 1.2)+
  geom_hline(yintercept=c(.95,0.75), color="black")+
  xlim(0,1.2)+
  geom_vline(xintercept = c(0.62,0.34))+labs(title="Base Model - EPA",
        x ="TP (mg/L)", y = "Probability of Occurence")

ggplot(Natural, aes(x=TP, col=Station_ID))+
  stat_ecdf(lwd=1.2)+
  geom_hline(yintercept=c(.95,0.75), color="black")+
  xlim(0,1.2)+labs(title="Model - Natural",
        x ="TP (mg/L)", y = "Probability of Occurence")

ggplot(Base_PS, aes(x=TP, col=Station_ID))+
  stat_ecdf(lwd=1.2)+
  geom_hline(yintercept=c(.95,0.75), color="black")+
  xlim(0,1.2)+labs(title="Base_PS",
        x ="TP (mg/L)", y = "Probability of Occurence")

TP concentrations and NP ratios visualized with cumulative distribution frequency

ggplot(BaseModel, aes(x = NP, col = Station_ID))+
   stat_ecdf(lwd = 1.2)+
  geom_hline(yintercept=c(.95,0.75), color="black")+
  xlim(0,10)+
  geom_vline(xintercept = c(5,3.1))+labs(title="Base Model - EPA",
        x ="NP Ratios", y = "Probability of Occurence")

ggplot(Natural, aes(x=NP, col=Station_ID))+
  stat_ecdf(lwd=1.2)+
  geom_hline(yintercept=c(.95,0.75), color="black")+
  xlim(0,10)+labs(title="Model - Natural",
        x ="NP Ratios", y = "Probability of Occurence")

LSPC188PERO.NoFrank<-subset(Base_NoFrank1,Station_ID=="LSPC188PERO")
dy.NoFrankTPLoad <- xts(x = LSPC188PERO.NoFrank$TPLoad, order.by = LSPC188PERO.NoFrank$Date)

LSPC188PERO.base<-subset(BaseModel1,Station_ID=="LSPC188PERO")
dy.TPLoad <- xts(x = LSPC188PERO.base$TPLoad, order.by = LSPC188PERO.base$Date)

LSPC188PERO.natural<-subset(Natural1,Station_ID=="LSPC188PERO")
dy.NatTPLoad <- xts(x = LSPC188PERO.natural$TPLoad, order.by = LSPC188PERO.natural$Date)

LSPC188PERO.Base_PS<-subset(Base_PS1,Station_ID=="LSPC188PERO")
dy.Base_PSTPLoad <- xts(x = LSPC188PERO.Base_PS$TPLoad, order.by = LSPC188PERO.Base_PS$Date)



LSPC188PERO.FrankPermit<-subset(Base_NoFrank1,Station_ID=="LSPC188PERO")
dy.FrankPermitTPLoad <- xts(x = LSPC188PERO.FrankPermit$TPLoad+174.5, order.by = LSPC188PERO.FrankPermit$Date)


#LSPC188PERO.Cont<-subset(Cont.w,Station_ID==c("usgs_03432350","USGS_0343233905"))
#dy.ContTPLoad <- xts(x = LSPC188PERO.Cont$TPLoad, order.by = LSPC188PERO.Cont$Date)

#LSPC188PERO.Grab<-subset(Grab.w,Station_ID==c("Frank_Up","Franklin_Site_1"))
#dy.GrabTPLoad <- xts(x = LSPC188PERO.Grab$TPLoad, order.by = LSPC188PERO.Grab$Date)


LSPC188PERO.TPLoad<-cbind(dy.TPLoad,dy.NatTPLoad,dy.Base_PSTPLoad,dy.NoFrankTPLoad,dy.FrankPermitTPLoad)


LSPC188PERO.TPLoad <- na.locf(LSPC188PERO.TPLoad)






p5 <- dygraph(LSPC188PERO.TPLoad, main="TP Load at LSPC188PERO - Below Franklin", ylab="TPLoad (lbs/day)") %>%
  dySeries("dy.TPLoad", drawPoints = TRUE, strokePattern="dashed", color="red") %>%
  dySeries("dy.NatTPLoad", drawPoints= T, color="blue") %>%
  dyOptions(labelsUTC = TRUE, fillGraph=TRUE, fillAlpha=0.1, drawGrid = FALSE, colors = RColorBrewer::brewer.pal(5, "Set1")) %>%
  dyRangeSelector() %>%
  dyCrosshair(direction = "vertical") %>%
  dyHighlight(highlightCircleSize = 5, highlightSeriesBackgroundAlpha = 0.2, hideOnMouseOut = FALSE)  %>%
  dyRoller(rollPeriod = 1)


p5

FrankPermitTPLoad = Permit limit 174 lbs/day added to No Franklin Model in order to represent TP Load in Harpeth with addition of Franklin Permit Limits

a<-subset(NPDES, Station_ID=="TN0020460")
b<-subset(NPDES, Station_ID=="TN0027278")
c<-subset(NPDES, Station_ID=="TN0028827")
d<-subset(NPDES, Station_ID=="TN0029718")
e<-subset(NPDES, Station_ID=="TN0057789")

f<-subset(NPDES, Station_ID=="TN0057827")
g<-subset(NPDES, Station_ID=="TN0057835")
h<-subset(NPDES, Station_ID=="TN0059790")
i<-subset(NPDES, Station_ID=="TN0060216")
j<-subset(NPDES, Station_ID=="TN0062332")

k<-subset(NPDES, Station_ID=="TN0064297")
l<-subset(NPDES, Station_ID=="TN0066958")
m<-subset(NPDES, Station_ID=="TN0067164")
n<-subset(NPDES, Station_ID=="TN0067873")
o<-subset(BaseModel, Station_ID=="HARPETH_1")


#c$TPLoad.csum<-c$TPLoad*0.00220462*NPDES$FLOW*86400

#MonthlyLoad <- function(x,y){
#  return(x*y*0.00220462*86400*30)
#}

#add_one_if <- function(vector, x_id, y_id, run_on_id){
#    if(vector[run_on_id]){
#        return(add_one(vector,x_id))}
#    else{
#        return(vector[x_id])
#    }
#}

#test$y <- apply(test, 1, add_one_if, x_id = 1, y_id = 2, run_on_id = 3)


#NPDES.list<-c("a","b","c","d","e","f","g","h","i","j","k","l","m","n")


#for (i in NPDES.list) {

#i$Date2=dplyr::lead(i$Date)

#i$Date.dif<-difftime(i$Date2,i$Date)


#}




#NPDES$TPLoad.csum <- ave(NPDES$TPLoad, NPDES$Station_ID, FUN=cumsum)

#NPDES$NP<-NPDES$TN/NPDES$TP

#NPDES1<-NPDES




#PS.TPload<-max(a$TPLoad.csum)+max(b$TPLoad.csum)+max(c$TPLoad.csum)+max(d$TPLoad.csum)+max(e$TPLoad.csum)+max(f$TPLoad.csum)+max(g$TPLoad.csum)+max(h$TPLoad.csum)+max(i$TPLoad.csum)+max(j$TPLoad.csum)+max(k$TPLoad.csum)+max(l$TPLoad.csum)+max(m$TPLoad.csum)+max(n$TPLoad.csum)

#PS.TPload #Cumulative TP Load between Jan 2011 and Dec 2019 of all combined Point Sources from NPDES permits
#max(o$TPLoad.csum) #Cumulative TP Load between Jan 2011 and Dec 2019 from HARPETH mouth (HARPETH_1)


#PS.TPload/max(o$TPLoad.csum)



#PS.TPload<-data.frame(c$Date,PS.TPload)

#colnames(PS.TPload)[1] <- "Date"
#colnames(PS.TPload)[2] <- "TPLoad.csum"

#p <- ggplot(data=PS.TPload, aes(x=Date, y=log(TPLoad.csum))) +
#  geom_line() + 
#  xlab("Year") +
#  geom_line(data=subset(BaseModel, Station_ID=="HARPETH_1"))

#p